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ABSTRACT 


Kink  instabilities  in  long  ion  layers  immersed  in  a  dense  back* 


ground  plasma  are  studied.  A  numerical  extension  of  the  analytic 
model  of  Lovelace  indicates  that  these  instabilities  will  occur  for  values 
of  the  self-magnetic  field  index  below  those  predicted  previously.  A 
quasineutral  hybrid  simulation  code  has  been  used  to  verify  these  lower 
thresholds.  The  simulations  also  show  that  the  end  of  exponential 
growth  occurs  due  to  a  nonlinear  shift  in  the  betatron  frequency  at  large 
amplitude,  producing  an  increase  in  layer  thickness  and  a  layer  which 
has  many  non-axis-encircling  ions. 
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I.  INTRODUCTION 

Field-reversing  ion  layers  have  been  considered  as  a  possible  means  of  confining 
and  heating  a  fusion  plasma.  It  has  also  been  proposed  that  such  ion  layers  could  pro¬ 
vide  a  stabilizing  influence  on  the  tilt  mode  in  the  spheromak1'2.  However,  ion  layers 
may  be  unstable  due  to  a  kink  mode  which  results  from  the  coupling  of  compressional 
Alfven  waves  with  the  betatron  motion  of  the  layer  particles. 

The  ion  layer  configuration  is  shown  in  Fig.  1.  An  ion  beam,  at  radius  R  and  of 
thickness  a,  flows  perpendicularly  to  an  external  magnetic  field,  Be—B£z.  The  beam 
current,  Jb9,  produces  a  self-magnetic  field,  £/,  which  may  be  large  enough  to  reverse 
the  total  magnetic  field  on  axis.  The  layer  is  immersed  in  a  background  plasma  such 
that  nb«  np,  where  nb  and  np  are  the  densities  of  the  beam  and  plasma,  respectively. 
The  background  plasma  is  assumed  to  be  sufficiently  dense  that  vj«  v} ,  where  vA  is 
the  background  plasma  Alfven  speed,  = B*/  (4v  np  Mp) 1/2,  and  v9  is  the  average  rota¬ 
tional  velocity  of  the  beam.  Mp  is  defined  to  be  the  mass  of  an  ion  in  the  background 
plasma.  A  conducting  wall  is  located  at  radius  r—rw.  There  is  no  variation  in  the  axial 
direction  (d/dz— 0).  It  is  this  case,  fc2—0,  which  is  expected  to  be  the  most  unstable 
for  a  long  layer.  This  has  been  indicated  by  Sudan  and  Rosenbluth3  through  the  appli¬ 
cation  of  an  energy  principle  derived  from  the  Vlasov  equation  to  obtain  sufficient 
conditions  for  stability.  Kink  modes  in  these  layers  correspond  to  perturbations  of  the 
form  T—  r* ,exp ( lm&—  iu t)  with  2.  This  paper  discusses  the  results  of  theory  and 
simulations  concerning  the  linear  and  nonlinear  behavior  of  these  modes.  The  term 
"kink”,  in  reference  to  these  modes,  arises  because  the  azimuthal  mode  number,  m,  is 
analogous  to  the  toroidal  mode  number  of  a  thin  ring,  or  "bicycle  tire",  configuration. 
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Section  II  presents  theoretical  results  for  instability  thresholds  and  linear  growth 
rates  as  determined  by  a  numerical  generalization  of  an  analysis  by  Lovelace4.  Section 
III  describes  the  model  for  our  hybrid  simulation  code  which  has  been  used  to  study 
ion  layer  kink  modes:  The  simulation  results  pertaining  to  the  linear  behavior  of  kink 
instabilities  are  presented  in  Sec.  IV.  Comparisons  are  made  with  the  results  of  Sec. 
n.  Section  V  describes  the  nonlinear  behavior  of  kink  modes  as  observed  in  the 
hybrid  simulations. 

II.  LINEAR  THEORY 

The  linear  stability  of  kink  modes  has  been  studied  by  Lovelace  in  the  approxi¬ 
mation  of  a  linear  beam5  6  and  as  a  long  layer4,  neglecting  axial  variation.  Similar 
instabilities  have  been  studied  by  Finn  and  Sudan7  in  the  bicycle  tire  limit.  Our  model 
follows  that  of  Lovelace4  in  which  kink  and  precessional  modes,  corresponding  to  per¬ 
turbations  in  radius  of  the  form  T-r«,exp(/mO— /«/),  were  studied.  In  this  paper  we 
consider  the  case  in  which  the  external  magnetic  field  is  uniform.  In  this  case  the  pre¬ 
cessional  mode  (m—  1)  is  neutrally  stable.  We  will  restrict  this  analysis  to  modes  with 
m>2  (i.e.,  kink  modes). 

The  analysis  of  Lovelace4  treated  the  background  plasma  ions  and  electrons  as 
fluids.  The  layer  was  assumed  to  be  composed  of  collisionless  ions  with  their  behavior 
described  by  the  particle  equations  of  motion.  The  layer  was  assumed  to  be  thin, 
(a/R)2«  1,  and  the  modes  were  assumed  to  be  of  low  frequency,  such  that 
Wn)2«l,  where  ft  is  the  layer  rotation  frequency.  This  assumption  allows  the 
forces  on  the  layer  particles  to  be  given  by  averages  over  the  layer  thickness.  The  per¬ 
turbation,  «r,  was  assumed  to  be  first  order,  (<r/a)«l,  as  were  the  perturbed  field 
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quantities,  SBZ/BZ«\,  8EJ Er«  1,  and  6E9/E9«  1.  The  equations  were  linearized; 
i.e.,  second  and  higher  order  terms  were  neglected.  The  following  equations  were 
obtained  by  Lovelace4  for  the  perturbation  e,: 


(<?/dthr--nH\-H,sUr  +  6Fr 


(la) 


^<8fl=>+«r<8£>-  a^<8£<> 

Mbc£l2  dr  z 

where  the  radial  averages  are  weighted  by  the  beam  current. 


(lb) 

(lc) 


<  (...)>  =  J  rdrJb(r)  (...)/ J  rdrJb(r). 

Mb  is  defined  to  be  the  mass  of  a  layer  ion,  v9  the  average  layer  rotation  velocity,  and 
r)s  is  the  self-magnetic  field  index,  which  is  approximately  equal  to  (to^/d)2.  The 
layer  betatron  frequency,  is  the  average  frequency  at  which  the  ions  oscillate  radi¬ 
ally  in  the  self-magnetic  field  of  the  layer.  Defining  (  to  be  the  field-reversal  factor, 

C=|(5/(r-0))/5/|, 

makes  r\s  of  order  (2/?/a){.  (Note  that  this  definition  of  {  is  not  the  same  as  used  by 
Lovelace4).  The  perturbed  field  quantities  are  determined  by  the  following  equations4: 


ICO 

c 

A* 

J-A2*r+4.| 

dr  r  r  | 

(f 

6E9  —  A  b ^(ABz+4>)  +  Afl 

dr  *  r 

where 


A,  3 


A*3 


1  +  (vVc)2 

(l+(vyc)2)2-(W«c/)2 


Jbj  _ (2 _ 

l  C  )  (l+CVcW-Ww*)2 

A*2S8*,-(4ir/c)cr/4 


(2a) 

(2b) 


(2c) 

(2d) 

(2e) 


0  =  (<or/  mv9)4irerJb/  c  , 

A  A  is  obtained  from  solving  the  differential  equation 


19  _  A  B  A  O  I  Iff |  9  A  )  A  a 

7d7("'-57iS'|  +  |“ — ~ JAB>"  T|  87a*|4B- 

- P7r _ ~«7 + 57 1  7  ) '  <3) 

Assuming  that  and  £#  vanish  at  a  conducting  wall  at  r—  rw,  Eq.  2b  provides  the 


boundary  condition  on  A9;: 


1  9  A  p 

A£r  9/  7 


imAh(rw) 


r—rm  A  a  (  rw) 


Note  that  A*  and  Aa  are  functions  of  radius  because  they  depend  on  vA  which,  in  a 
strong  ion  layer,  will  normally  depend  on  r. 

Lovelace4  used  Eqs.  1,2,3,  and  4  to  solve  for  the  instability  thresholds  and 
growth  rates.  In  order  to  do  this  analytically,  it  was  assumed  that  vj=constant  and 
that  the  layer  had  sharp  boundaries.  Thresholds  were  obtained  for  two  cases:  (1) 
at2«(vA/ a)2  and  (2)  a>2<  ( vA/ a)2  with  a>2«y2,  where  <ur  and  y  are  the  real  and 
imaginary  parts  of  the  frequency.  In  both  cases  17  *<  n^—\  was  found  to  be  sufficient 
for  stability.  The  assumption  v^=constant  implies  that  the  plasma  density  is  smaller 
inside  than  outside  the  layer.  This  assumption  does  not  correspond  to  the  situation  in 
real  three-dimensional  systems  when  the  field-reversal  factor  is  close  to  or  greater  than 
unity.  This  is  because  in  real  systems  the  density  will  not  have  large  variations  along  a 
given  closed  field  line.  The  comment  is  made  by  Lovelace4  that  it  appears  feasible  to 


solve  Eq.  3  numerically  in  order  to  remove  this  difficulty. 


We  have  chosen  to  extend  the  analysis  of  Lovelace4  in  order  to  eliminate  the 
assumption  of  constant  Alfven  speed,  as  well  as  the  frequency  restrictions  m2<(vA/a)2 
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and  ot}<y2>  by  performing  a  numerical  solution  of  the  system  of  Eqs.  1,  2,  3,  and  4. 
The  procedure  used  for  this  numerical  solution  is  described  in  Appendix  A.l.  Rather 
than  a  sharp  boundary  layer  equilibrium  distribution,  an  exponential  rigid  rotor  equili¬ 
brium8  is  treated.  The  beam  current  density  for  this  distribution  is 


J w (r)  —  eriofl  rsech2 


r 2-r? 


and  the  magnetic  field  profile  is 


.  -2Mbv}c  |  r^-r2) 

Bt(r)  -  — tanh  — r-M  (6) 

eSlri  (  r$  j 

In  Eqs.  5  and  6,  v,  is  the  thermal  velocity  of  the  rigid  rotor  distribution,  rx  and  r0  are 
adjustable  parameters,  and  n0  is  a  normalization  factor  for  the  density.  The  equili¬ 


brium  background  plasma  density  is  considered  to  be  uniform. 


Our  numerical  results  show  significant  differences  from  the  analytic  results  of 
Lovelace4.  In  particular  we  find  the  instability  threshold  for  a  given  mode  to  occur  at 
somewhat  lower  values  of  t)s  than  m2— 1.  Additionally,  we  have  found  that  the  insta¬ 
bility  thresholds  are  influenced  by  the  Alfven  transit  time  across  the  layer  thickness. 
The  approximate  thresholds  are  summarized  in  Table  1  for  different  values  of  the 
Alfven  transit  time.  One  can  see  from  Table  1  that  for  a  fixed  layer  thickness,  a,  an 
increased  background  density  provides  a  stabilizing  effect.  The  mode  with  the  lowest 
instability  threshold  in  i),  is  the  m—2  mode.  The  low  thresholds  here  imply  that  a 
thick  layer,  a~R,  is  a  necessary  condition  for  an  ion  layer  to  be  stable  near  field- 
reversal  ({=1)  (this  condition  can  only  be  called  a  necessary  one  since  this  analysis 
requires  a«r  to  determine  a  growth  rate).  These  results  were  not  available  under 
the  simplifying  assumptions  used  by  Lovelace4  because  near  the  thresholds  «,>y  and 
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The  real  and  imaginary  parts  of  the  frequency  of  modes  for  a  layer  with  (fl/a)-5 
are  shown  as  a  function  of  field-reverse  factor,  £,  in  Figs.  2  and  3.  The  maximum 
growth  rate  for  a  given  mode  occurs  when  Lovelace4  found  the 

growth  rates  to  be  linearly  proportional  to  the  Alfven  speed.  We  find  that  the  critical 
factor  is  not  simply  the  Alfven  speed,  but  the  Alfven  transit  time  across  the  layer 
thickness.  The  growth  rate  of  the  w— 2  mode  for  layers  of  varying  background  density 
are  shown  in  Fig.  4  for  a  layer  with  R/as »5.  As  the  Alfven  transit  time  becomes  very 
large  and  a>2»  vj/  a2,  then  the  mode  will  be  stabilized  because  the  layer  will  no 
longer  behave  rigidly. 

For  cases  in  which  the  assumptions  of  Lovelace4  are  valid,  our  growth  rates 
reduce  to  those  obtained  analytically. 

III.  SIMULATION  MODEL 

The  simulation  code  that  we  have  used  to  study  ion  layer  kink  instabilities  is  a 
two-dimensional,  fully  nonlinear,  quasineutral  hybrid  code.  Such  codes  have  been 
used  in  one-dimensional  problems  by  Byers  et  al.9  and  in  two-dimensional  studies  of 
theta  pinch  implosions  by  Hewett10.  Because  ion  layers  are  characterized  by  ion  Lar- 
mor  radii  comparable  to  the  system  size,  it  is  necessary  to  treat  the  ions  as  particles. 
Since  the  time  scales  of  kink  modes  are  much  longer  than  the  ion-cyclotron  period,  it 
is  unnecessary  to  follow  the  details  of  the  electron  dynamics.  Hence,  the  electrons  are 
considered  to  be  an  inertialess  fluid.  In  addition,  electromagnetic  radiation  effects  are 
not  important  on  these  time  scales,  allowing  us  to  use  the  Darwin  version  of 
Ampere's  law  (i.e.,  the  transverse  displacement  current  is  neglected). 

The  Darwin  version  of  Ampere's  law  is  combined  with  the  inertialess  version  of 


8 


the  electron  momentum  equation  and  the  assumption  of  quasineutrality  to  determine 
the  electric  field.  The  expression  for  the  electric  field  (as  derived  in  Ref.  11,  or  simi¬ 
larly  in  Ref.  9)  is 

£--~—(Vxfl)xg--±-7xg-  — V(rt,Tf).  (7) 

4 itn,e  it/ec  n,e 

In  Eq.  7  the  ion  density,  nh  and  ion  current  density,  7„  are  determined  from  the  ion 
particles  by  linear  weighting  (particle-in-cell)  from  the  grid.  The  magnetic  field  is 
advanced  in  time  by  Faraday’s  law  and  the  ions  are  moved  by  the  equations  of 
motion.  A  more  detailed  description  of  this  code  may  be  found  in  Ref.  11. 

The  equations  are  solved  in  the  r— 9  plane  (see  Fig.  1).  However,  for  numerical 
reasons,  cartesian  coordinates  are  used  in  the  simulation.  No  axial  variation  is  allowed 
G.e.,  d/dz— 0)  and  Bx,  By,  £.,  Jz,  and  v;  are  all  set  equal  to  zero.  Doubly  periodic 
boundaries  have  been  used  in  the  simulation  results  discussed  here.  A  conducting 
wall  at  radius  rw  has  a  weakly  stabilizing  effect  on  ion  layer  kink  modes.  The  results 
of  Lovelace4  indicate  that  the  effect  of  the  wall  will  be  negligible  unless 
( rw—R)/R  « 1;  we  have  confirmed  this  result  by  the  numerical  calculations  described 
in  Sec.  U. 


Our  simulations  begin  with  the  loading  of  a  cold  uniform  background  plasma  of 
fluid  electrons  and  particle  ions  with  a  uniform  density  np.  Then  a  Vlasov  equilibrium 
(i.e.,  d/o/d  f—0)  ion  layer  distribution,  corresponding  to  an  exponential  rigid  rotor, 


/(r,  vr,  ve)  -  -~sech2 
wv,2 


.-m)2]] 

2 


(8) 


is  added.  This  distribution  produces  the  ion  current  and  magnetic  field  profiles 


described  in  Sec.  n.  Fluid  electrons  are  present  in  the  layer,  in  the  same  initial  profile, 
to  provide  charge  neutralization.  However,  all  current  is  initially  carried  by  the  ions. 


The  ion  layer  is  assumed  to  be  tenuous  relative  to  the  background  plasma,  so  that 
nb«np  and  vj}«Tg,  as  in  Sec.  II.  The  layer  ions  are  represented  by  40,000  parti¬ 
cles.  40,000  particles  are  also  used  to  represent  the  background  plasma  ions.  Finite 
differencing  is  performed  on  a  100  by  100  grid  and  the  time  step  is  <wc,A  r— 0.1. 

IV.  SIMULATION  RESULTS:  LINEAR  BEHAVIOR 

Simulations  have  been  performed  using  the  code  described  in  Sec.  Ill  to  study 
the  linear  and  nonlinear  behavior  of  ion  layer  kink  instabilities.  Figure  5  shows  the 
initial  particle  positions  for  an  exponential  rigid  rotor  ion  layer;  the  background  plasma 
is  also  represented  by  particles,  but  these  are  not  shown  on  this  plot.  The  field- 
reversal  factor  for  this  case  is  (—1.1  ( B .  slightly  reversed  on  axis).  The  rigid  rotor  dis¬ 
tribution  of  Eq.  6  has  parameters  set  to  ri— 2.5  and  r0— 1.35,  corresponding  to  R/a=6. 
The  self-magnetic  field  index  for  this  case  is  i)s- 5.9  and  the  normalized  inverse 
Alfven  transit  time  is  v^/w^a^O.I.  Table  I  indicates  that  this  Vlasov  equilibrium 
should  be  unstable  to  both  the  m— 2  and  m— 3  modes.  Figure  6  shows  the  layer  at 
time  /— 50«dl,  when  an  m— 3  mode  has  grown  to  large  amplitude.  The  m— 3  character 
is  clearly  evident  in  Fig.  7,  which  shows  the  azimuthal  electric  field  contours.  At  a 
later  time,  t— 100a>d',  an  m— 2  mode  clearly  dominates,  as  can  be  seen  in  Figs.  8  and 
9.  At  a  much  later  time,  r—  SOOw^1,  the  layer  has  become  very  thick  and  the  earlier 
mode  structure  has  disappeared,  as  seen  in  Fig.  10. 

The  numerical  calculation  described  in  Sec.  II  predicts  the  linear  growth  rate  for 
the  m— 2  mode  to  be  y/a>d— 0.065,  and  for  the  /n— 3  mode  to  be  y/«d— 0.112. 
Growth  rates  were  obtained  from  the  simulation  by  measuring  the  perturbation  in  the 
average  radial  particle  position  as  a  function  of  angle,  8r(0,r)— (<  r(0, /)>—/?).  Br(B) 


10 


is  then  decomposed  into  its  Fourier  components.  (8r(r))2  for  the  m— 2  and  m— 3 
modes  is  shown  in  Figs.  11  and  12.  The  growth  rate  for  the  m— 2  mode  is  found  to  be 
y/<DC(— 0.065,  and  the  growth  rate  for  the  m— 3  mode  is  found  to  be  y/<yC(— 0.101,  very 
close  to  the  expected  values.  The  m—\  mode  and  modes  with  m> 4  were  found  to  be 
stable,  as  predicted.  In  Fig.  2  the  growth  rates  measured  from  simulations  are  plotted 
as  a  function  of  field-reversal  factor,  with  the  results  from  Sec.  II  for  modes  m—2,  3, 
4,  and  5.  Figure  4  shows  the  results  of  simulation  and  theory,  giving  growth  rates  as  a 
function  of  the  inverse  Alfven  transit  time.  In  both  cases  we  find  good  agreement  for 
linear  growth  rates.  The  simulations  have  also  confirmed  the  predictions  that  the  ins¬ 
tability  thresholds  should  be  lower  in  terms  of  tj5  than  the  previous  analytic  results4. 

In  order  to  obtain  the  azimuthal  mode  components  of  (8r)2  (shown  in  Figs.  11 
and  12)  used  to  determine  the  growth  rates,  the  values  of  8  r(0)  from  0— 0  to  0—  2ir 
were  used.  Because  this  information  always  contains  an  integer  number  of 
wavelengths  for  an  unstable  mode,  the  real  frequency  of  oscillation  cannot  be  obtained 
from  the  growth  rate  diagnostics.  Real  frequencies  were  determined  by  measuring  the 
frequency  of  the  radial  perturbation  at  a  fixed  angle.  Because  the  real  frequency  of 
these  modes  is  comparable  to  the  growth  rate,  (Sr)2  may  grow  several  decades  in  one 
period  of  oscillation.  This  made  the  measurement  of  real  frequencies  difficult  except 
for  cases  in  which  one  mode  dominated  for  a  long  period  of  time.  For  the  limited 
number  of  cases  in  which  real  frequencies  could  be  evaluated,  they  were  found  to 
agree  with  the  predictions  of  Sec.  II  to  within  ten  percent. 

The  physical  mechanism  for  ion  layer  kink  instabilities  is  the  resonant  interaction 
of  the  betatron  motion  of  layer  particles  with  compressions!  Alfven  waves  in  the  back¬ 
ground  plasma,  which  occurs  when  (o=a>^—mCl.  The  maximum  growth  rates  are 
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expected  to  occur  when  the  largest  number  of  particles  are  experiencing  resonant  beta¬ 
tron  motion.  This  should  occur  when  since  <u2«fi2.  This  has  been 

observed  in  both  our  simulations  and  in  the  numerical  calculations  described  in  Sec. 
II,  as  the  peaks  in  Fig.  2  correspond  to  points  where  The  growth  of  the 

m— 3  mode  in  the  preceding  example  may  be  traced  to  the  fact  that  most  particles  in 
this  layer  have  betatron  motion  characterized  by  3.  An  example  of  a  particle 

orbit  from  f—0  to  /— 60o>^'  is  shown  in  Fig.  13,  with  the  m— 3  character  clearly 
present.  One  reason  that  thresholds  may  occur  at  lower  values  of  the  self-magnetic 
field  index  than  7)s—  m2—  1  is  because  even  when  the  mean  betatron  frequency  of  the 
layer  is  too  small  to  produce  a  resonant  interaction,  the  thermal  spread  of  the  layer 
ions  can  still  allow  a  substantial  number  of  particles  to  resonate  at  a  higher  betatron 
frequency  and  drive  the  instability. 

V.  SIMULATION  RESULTS:  NONLINEAR  BEHAVIOR 

Nonlinear  effects  which  halt  the  growth  of  ion  layer  kink  instabilities  have  been 
indentified.  Figures  11  and  12  show  the  end  of  exponential  growth,  or  saturation  of 
the  w— 2  and  m— 3  instabilities  for  the  unstable  configuration  of  Sec.  IV.  The  common 
observation  at  saturation  has  been  a  nonlinear  shift  in  the  layer  betatron  frequency 
which  diminishes  the  resonant  interaction  of  layer  particles  required  to  produce  kink 
instabilities. 

During  the  growth  of  kink  instabilities,  a  large  fraction  of  the  directed  beam 
energy  is  lost  to  thermal  energy.  This  can  be  seen  in  Figs.  14  and  15,  which  show  the 
directed  and  thermal  energies  of  the  layer  as  a  function  of  time.  The  thermal  layer 
pressure  increases  and  is  no  longer  in  equilibrium  with  the  magnetic  forces  compress* 
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ing  the  layer.  This  pressure  imbalance  causes  the  layer  to  expand.  The  change  in 
layer  thickness  is  clearly  visible  in  Figs.  6  and  8,  after  the  m— 2  and  m— 3  instabilities 
have  grown  to  large  amplitude.  A  measure  of  the  expansion  with  time  is  shown  in 
Fig.  16,  which  plots  <(r—  <r(0)>)2>  as  a  function  of  time.  The  significant  point  is 
that  the  layer  betatron  frequency  is  roughly  proportional  to  R/a.  Therefore,  as  the 
layer  increases  its  thickness,  a,  the  betatron  frequency  and  consequently  the  self- 
magnetic  field  index  is  reduced.  When  the  layer  becomes  sufficiently  thick  so  that  the 
self-magnetic  field  index  drops  below  the  instability  threshold  for  a  given  mode,  that 
mode  ceases  to  grow  exponentially,  although  nonlinear  effects,  producing  layer  heat¬ 
ing,  have  been  observed  to  persist  somewhat  beyond  that  point.  An  additional  factor 
is  that  the  layer  expansion  increases  the  Alfven  transit  time  across  the  layer,  which 
tends  to  stabilize  the  instability.  This  effect,  however,  is  small  compared  to  that  pro¬ 
duced  by  the  shift  in  betatron  frequency. 

The  effects  of  the  nonlinear  shift  in  betatron  frequency  can  be  seen  in  the  exam¬ 
ple  of  Sec.  IV.  The  m— 3  mode  has  the  largest  growth  rate  and  is  the  dominant  mode 
at  early  times.  However,  at  n-560*^1  the  m— 3  mode  stops  growing.  At  this  point  the 
self-magnetic  field  index  is  4.35  with  0.074.  It  can  be  seen  from  Table  1 

that  these  values  are  very  close  to  the  linear  threshold  of  the  m— 3  mode.  The  m— 2 
mode  is  still  linearly  unstable  and  continues  to  grow  until  it  eventually  exceeds  the 
amplitude  of  the  saturated  m—3  mode,  and  then  becomes  the  dominant  mode,  as  seen 
in  Fig.  8.  At  r— the  m— 2  mode  stops  growing.  At  this  point  vj,— 2.1 1  and 
vA/a<aci~ 0.036.  The  Table  1  values  can  only  be  expected  to  give  an  approximate 
result  for  this  mode,  since  the  layer  has  become  thick  enough  so  that  the  self- 
magnetic  field  index  is  only  an  approximate  concept  as  defined  by  Eq.  lc.  Neverthe- 


less,  the  m— 2  mode  stops  growing  at  a  time  when  when  t),  is  very  near  the  threshold 
predicted  by  Table  1.  Similar  cascades  have  been  observed  for  higher  azimuthal  mode 
numbers.  For  higher  azimuthal  mode  numbers  the  mode  with  the  largest  mode 
number  will  saturate  first,  followed  by  the  successive  saturation  of  the  lower  mode 
numbers  beginning  with  the  largest  remaining  unstable  mode  and  ending  with  the 
saturation  of  the  m— 2  mode.  It  is  rarely  possible,  however,  to  see  clearly  more  than 
two  modes  because  the  layer  expands  rapidly  enough  that  the  betatron  frequency  will 
drop  below  the  thresholds  of  modes  before  they  can  grow  to  a  sufficiently  large  ampli¬ 
tude  to  dominate.  Nevertheless,  the  final  state  is  similar,  regardless  of  the  number  of 
unstable  modes  present  in  the  initial  equilibrium.  This  has  been  confirmed  in  simula¬ 
tions  of  thinner  layers,  which  have  had  modes  up  to  m— 5  unstable,  yet  result  in  virtu¬ 
ally  the  same  final  state  as  the  previously  discussed  example. 

After  all  of  the  initial  instabilities  have  saturated,  the  resultant  thick  layer  appears 
to  be  stable.  No  further  instabilities  have  been  observed  for  times  as  long  as 
500a)"1.  However,  the  final  state  has  a  very  high  noise  level  and  we  cannot  exclude 
the  possibility  that  instabilities  may  exist  with  smaller  growth  rates  than  those  of  the 
initial  kink  modes.  The  layer  is  characterized  by  being  thick,  but  still  field-reversed. 
This  can  be  seen  by  the  initial  and  final  magnetic  field  profiles  shown  in  Fig.  17.  For 
stronger  layers,  with  field-reversal  factors  as  large  as  {—1.8,  the  observed  behavior  has 
been  qualitatively  the  same,  in  that  the  final  state  is  much  thicker,  relatively  stable, 
and  field-reversal  is  maintained.  The  instabilities  do  not  result  in  the  loss  of  the  layer, 
but  the  final  state  has  substantial  electron  currents  and  many  non-axis-encircling  ions. 
The  nature  of  the  change  in  orbits  can  be  seen  in  Fig.  18.  The  large  axis-encircling 
part  of  the  orbit  (see  Fig.  13)  corresponds  to  the  initial  m—3  betatron  motion,  which 
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drives  the  initial  instability.  As  the  layer  broadens,  the  motion  changes  to  m— 2  beta¬ 
tron  motion  and  finally  the  orbit  becomes  non-axis-encircling  when  the  layer  thickness 
becomes  comparable  to  its  major  radius. 

As  a  check  on  numerical  accuracy,  total  energy  has  been  monitored  and  has  been 
found  to  be  conserved  to  within  two  percent  of  the  variation  of  potential  and  kinetic 
energies.  Additional  results  of  ion  layer  kink  instability  simulations  are  presented  in 
Appendix  A.  2. 

VI.  SUMMARY 

The  analysis  of  Lovelace4  has  been  generalized  numerically.  The  numerical 
results  show  that  the  thresholds  are  at  lower  values  of  the  self-magnetic  field  index 
than  predicted  by  earlier  analytic  results.  Simulations  have  confirmed  these  lower 
values  of  the  instability  thresholds.  Instabilities  have  been  found  to  grow  on  magne¬ 
tohydrodynamic  time  scales;  i.e.  y^vja.  Nonlinear  effects  have  been  found  to 
reduce  the  betatron  frequency  of  the  layer  so  that  ion  betatron  motion  is  no  longer 
resonant  with  the  wave,  halting  exponential  instability  growth. 

Although  this  study  is  only  two-dimensional,  the  same  type  of  behavior  might  be 
expected  for  the  fully  three-dimensional  case,  since  Finn  and  Sudan7  have  obtained 
similar  results  for  the  linear  stability  of  ion  layers  in  the  bicycle  tire  limit.  It  is  possible 
that  the  kink  instability  could  be  stabilized  by  the  addition  of  a  toroidal  component  of 
the  magnetic  field12.  Without  a  toroidal  field  it  appears  that  the  use  of  ion  layers,  car¬ 
rying  a  substantial  part  of  the  current  needed  to  produce  field-reversal,  is  possible  only 
for  thick  layers  with  large  numbers  of  non-axis-encircling  particles. 
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Instability  Thresholds 
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Q. 

Vs 

m- 2 

m— 3 

m=4 

m=5 

.20 

0.4 

2.2 

5.0 

8.5 

.10 

1.0 

4.2 

8.0 

12.0 

.05 

1.8 

5.0 

9.5 

15.0 

.025 

2.0 

6.0 

11.0 

18.0 

m2-! 

3.0 

8.6  1 
_ 1 

15.0 

24.0 

TABLE  1.  Approximate  ion  layer  kink  instability  threshold  values  of  the  self-magnetic 

field  index,  tj,,  at  different  values  of  (vA/uela).  The  thresholds  were  determined  by 

% 

the  numerical  calculations  discussed  in  Sec.  II.  The  threshold  values  of  tjj  as  deter¬ 
mined  analytically  by  Lovelace,  1,  are  also  given. 
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FIG.  2.  Linear  growth  rates,  y,  for  ion  layer  kink  instabilities  as  a  function  of  field- 
reversal  factor,  The  solid  curves  are  the  results  of  the  numerical 

solution  of  Sec.  II.  The  points  are  the  results  of  hybrid  simulations.  R/a—S  and 
vA/(awci)~ 0.1. 
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FIG.  14.  Directed  beam  energy,  where  the  summation  is  over  all  beam 

i 

particles,  as  a  Function  of  time,  showing  the  loss  of  directed  energy  during  instability 
growth. 
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FIG.  15.  Beam  thermal  energy,  V2+v*,l»  where  the  summation  is 

i  2  ' 

over  all  particles,  as  a  function  of  time,  showing  the  beam  heating  during  instability 
growth. 


3 


FIG.  16.  <(r—<r(0)>)2>  as  a  function  of  time,  showing  the  layer  expansion.  The 
expansion  from  f-0  to  <-6 is  due  to  the  m- 3  instability.  At  later  times,  from 
r— to  /— 140w“',  the  expansion  is  due  to  the  m— 2  insubility. 
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FIG.  17.  Magnetic  Held  profiles,  Bf{r)  averaged  over  0,  at  r— 0  (dotted  line)  and  at 
r*300ttf^],  after  instabilities  have  saturated  (solid  line). 
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A.l  Numerical  Calculation  of  Linear  Growth  Rates  for  Ion 

Layer  Kink  Instabilities 


This  appendix  describes  the  procedure  used  to  calculate  numerically  the  complex 
frequencies  from  the  equations  of  Lovelace1.  Assuming  perturbations  of  the  form 
7,(0,  r)— re rexp(/m0—  m  t)  Eq.  la  can  be  written  as 

(m2ft2— 2mnw+ci>2)€,>-—  02(l+7ji)e,  +  8F,.  (A.l) 

The  radial  perturbation  is  rigid,  i.e.  «,(0,f).  With  this  assumption,  Eq.  la  then 

becomes 


f(<n)  •  —  m2llJ+2mftoi- <u2—  ft2(l+ijj)  +  — -87,  “0  (A.2) 

€r 

The  object  is  to  find  the  complex  roots,  o»,  of  Eq.  A.2.  This  process  would  be  straight* 
forward,  except  for  the  fact  that  8Fr  is  a  complicated  function  of  a*,  involving  integra¬ 
tion  over  the  radial  thickness  of  the  layer  and  the  solution  of  a  second  order 
differential  equation. 


If  a  guess  is  made  for  the  complex  frequency,  the  following  method  will  deter¬ 
mine  the  corresponding  value  of  (1  /tr)8Fr  With  the  known  rigid  rotor  profiles  for 
the  current  density  and  magnetic  field,  one  can  obtain  A a(r),  A b(r),  and  4'(r)  from 


.  | T  l+uyc)1 

"(c)  (l+(»,/d!)!-W«.„)1 

A  [  v*  _ its. _ 

M  c  )  (l+(vVc)2)2-(a./w„)2 


(A.3a) 

(A.3b) 
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4>'=-&—  JUJ-UnJtJc 

*r  |mv9j 

It  is  then  possible  to  solve  for  A  B'z  from  the  second  order  differential  equation 

if  L,  f A  B'.j  +  4-  -  -|f  A  »U', 

r  dr  dr  c2  r  r  dr 

m2Aa  (omve\^  1  0  ^  ,  d  I  imAb4>' 

“  — 5 - 5 — <t> - — lAa0 )  +  - . 

r2  c2r  )  r  dr  dr|  r 

At  the  conducting  wall  at  r—  rw  we  have  the  boundary  condition 


(A.  3c) 


1  Bad1 

Ss7a74i- 


imAb(rw ) 
rw  A  n  ( rw ) 


and  at  r— 0,  A2TZ(0)— 0.  This  is  a  second  order  two-point  boundary  value  problem 
and  is  solved  by  finite  difference  methods.  Once  A  B':  is  determined  it  is  possible  to 
obtain  85.,  8  En  and  8£fl  from 


(A.  6a) 


J2L**L - A  a  —(A 5\-h*')  “  A  *  |-A5\.+  £- 


C  «r 


i2LMk._-A6-i2L(A5'z+0')+Aa  |-A5'.+£- 

c  e,  r  if  ‘  r 


(A.6b) 


(A.6c) 


The  quantities  ij*,  <85.>,  <  8  £,  > ,  and  <8£#>,  which  are  used  in  8£,,  are 
evaluated  by  numerical  integration  over  the  layer  using  expression  for  the  average, 


<(...)>  =  J* rdrJb(r) (...)/ J  rdrJb(r)  . 


Then,  using 


—  <«5r>  +  -J--L<8£r>  - 

«r  Mbc*r  Mb  €r  Mbm€r 


<8£#> ,  (A.  /) 


all  of  the  terms  in  Eq.  A.2  are  now  known  for  the  given  frequency  guess. 


The  above  procedure  is  repeated  for  each  guess  of  u  as  a  function  call  in  a 
Muller’s  method2  subroutine.  Muller’s  method  will  provide  guesses  for  the  complex 
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roots  of  /(«)  in  Eq.  A.2  until  the  given  convergence  criteria  are  met.  The  Muller’s 
method  subroutine  used  in  this  code  is  based  on  a  subroutine  written  by  Au-Yeung 
and  Friedman3. 
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A.2  Further  Results  of  Ion  Layer  Kink  Instability  Simulations 


This  appendix  presents  results  from  two  runs  which  demonstrate  the  effects  of 
increasing  field-reversal  and  decreasing  layer  thickness  on  ion  layer  kink  instabilities. 

1.  Decreased  Layer  Thickness 

In  this  simulation  we  use  the  same  background  density  and  field-reversal  as  used 
in  the  run  described  in  detail  in  the  main  section,  i.e.,  »,,/««„- 0.1  and  {—1.1.  How¬ 
ever,  we  reduce  the  initial  layer  thickness  so  that  10  rather  than  6.  This  has  the 
effect  of  increasing  the  self-magnetic  field  index  from  tjs— 5.9  to  tj 1 3.0.  This  value 
of  t)s  is  slightly  larger  than  the  threshold  for  the  m— 5  mode.  Therefore,  we  expect 
modes  2,  3,  4,  and  5  to  be  unstable.  Our  simulation  results  confirm  this.  The  growth 
of  the  m—5  mode  is  relatively  small  and  is  not  evident  in  the  snapshot  diagnostics. 
Figure  1  shows  the  initial  particle  positions  for  this  thin  layer.  In  Fig.  2,  at  t-30oj~\ 
an  m— 4  instability  is  apparent,  as  expected,  since  its  growth  rate  is  the  largest  of  the 
unstable  modes.  This  mode  ceases  to  grow  at  r— The  m— 3  mode,  having  a 
smaller  growth  rate  than  the  m-4  mode,  continues  to  grow  and  its  amplitude  eventu¬ 
ally  exceeds  the  amplitude  of  the  m— 4  mode.  This  can  be  seen  in  Fig.  3  at  r— TOw^1. 
At  a  considerably  later  time,  J— the  m-2  mode  dominates,  as  is  clear  from 
Fig.  4.  The  time  evolution  of  (Sr)2  is  plotted  in  Figs.  5-8  for  these  modes.  Figure  9 
shows  the  change  in  <  (r-<  r(B)>)2>  with  time.  The  increase  in  layer  thickness  due 
to  each  of  the  modes,  m-2,  3,  and  4,  can  be  seen  in  this  figure.  Although  this  run  is 
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not  as  long  as  the  one  described  in  the  main  section,  note  that  at  r— 220a>"1,  after  all 
modes  have  saturated,  the  magnitude  of  <(/■— <r(0)>)2>  at  a  comparable  time  is 
approximately  the  same,  indicating  that  the  only  effect  of  making  a  layer  thinner  is 
that  it  must  proceed  through  the  saturation  of  more  instabilities  before  it  reaches  the 
eventual  thick  layer  state. 

2.  Increased  Field-Reversal 

This  run  has  the  same  background  plasma  density  as  the  run  discussed  in  the 
main  section.  The  layer  thickness  is  slightly  larger,  Vt/a-5  rather  than  6.  The  pri¬ 
mary  difference  is  that  the  degree  of  field-reversal  has  been  increased  from  {—1.1  to 
(—1.38.  This  has  the  effect  of  increasing  the  self-magnetic  field  index  from  t)5— 5.9  to 
t|j— 12.3.  This  value  is  sufficiently  large  so  that  the  m— 4  mode  will  be  unstable  in 
addition  to  the  m— 2  and  m— 3  modes.  Figure  10  shows  the  initial  particle  positions 
for  this  layer.  Figure  11  shows  the  layer  at  t— SOu^1,  after  and  m— 4  instability  has 
grown  to  large  amplitude.  In  Fig.  12,  at  t~7(kj~l,  the  m— 4  mode  has  saturated  and 
the  m— 3  mode  now  dominates.  Finally,  in  Fig.  13,  at  r— 220<u~1,  after  the  exponential 
growth  of  all  modes  has  stopped,  some  m-2  structure  still  remains  due  to  the  growth 
of  the  last  unstable  mode,  the  m— 2  mode.  This  run  ended  at  this  point;  however, 
from  the  results  of  other  similar  runs,  such  as  the  one  described  in  the  main  section, 
we  would  expect  this  structure  to  eventually  dissappear.  Figure  14  shows  the  change 
in  the  magnetic  field  profile.  Note  that  field-reversal  is  maintained. 
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FIG.A3  Particle  positions  at  The  layer  has  thickened  and  now  an  m— 3  ins¬ 

tability  dominates. 
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FIG.  A6  (8  r)3  as  a  function  of  time  for  the  3  mode. 
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FIG.  a 7  (8r)2  as  a  function  of  time  for  the  m- 4  mode. 


FIG.  a9  <  (r— <  r(0)>)2>  as  a  function  of  time. 
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FIG. All  Particle  positions  at  /-50o)~'.  An  m-4  mode  has  grown  to  large  amplitude. 


FIG.A12  Particle  positions  at  An  m— 3  mode  now  dominates. 
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FIG.A13  Particle  positions  at  r— 220tu^1.  The  structure  of  the  /n- 2  mode,  the  last 
mode  to  saturate,  is  still  apparent. 


